# Cleans the results from fitting PLANT elec gen model
  library(pacman)
  p_load(
    here, fst, data.table, purrr, janitor, ggplot2, lubridate, stringr,
    scales
  )

  # Regional load data
  nerc_load_dt = 
    read.fst(
      path = here("Data/electricity-generation/clean-load/nerc-load-dt-oge.fst"),
      as.data.table = TRUE
    )[datetime_utc %between% 
      c(ymd_hms('2019-01-01 10:00:00'),ymd_hms('2019-12-31 23:00:00'))
    ][,
      ic := factor(fcase(
        nerc_adj %in% c('CAL','WECC'), 'West',
        nerc_adj == 'TRE', 'Texas',
        default = 'East'
      ), levels = c('West','Texas','East'))
    ]
  # Plant info table
  plant_info_dt = 
    read.fst(
      path = here("Data/electricity-generation/plant-info-dt-oge.fst"),
      as.data.table = TRUE
    )[year == 2019,.SD[which.min(year)], by= plant_id_eia] |>
    setkey(plant_id_eia, year)
  # Reading the electricity generation files (with fitted values)
  # Need to group by ic, nerc, fuel category and time 
  summarize_plant_gen_dt_nerc = function(nerc_adj_in, plant_info_dt){
    # Getting file paths for all plants in region
    plant_gen_fp = list.files(
      here('Data/electricity-generation/plant-model-fit/plant-gen-dt'),
      full.names = TRUE
    )
    nerc_plant_gen_fp = 
      plant_gen_fp[
        str_extract(plant_gen_fp, '(?<=-)\\d*(?=\\.fst)') |> as.integer()
        %in% plant_info_dt[nerc_adj == nerc_adj_in]$plant_id_eia
      ]
    # Summarizing generation by time, ic, nerc, fuel cat
    elec_gen_dt = 
      map_dfr(
        nerc_plant_gen_fp, 
        \(fp_in){
          read.fst(
            path = fp_in,
            as.data.table = TRUE
          ) |>
          merge(
            plant_info_dt, 
            by = 'plant_id_eia'
          )
        }
      )[,.(
        tot_net_generation_mwh = sum(net_generation_mwh),
        tot_fit_net_generation_mwh = sum(fit_net_generation_mwh)), 
        keyby = .(
          datetime_utc,
          ic = factor(fcase(
              nerc_adj %in% c('CAL','WECC'), 'West',
              nerc_adj == 'TRE', 'Texas',
              default = 'East'
            ), levels = c('West','Texas','East')),
          nerc_adj, 
          fuel_category
        )
      ]
        
    return(elec_gen_dt)
  }

  elec_gen_dt = 
    map_dfr(
      c('CAL','MRO','NPCC','RFC','SERC','TRE','WECC'),
      summarize_plant_gen_dt_nerc, 
      plant_info_dt = plant_info_dt
    )


# Summarizing production by region
nerc_gen_dt = 
  elec_gen_dt[!(fuel_category %in% c('solar','wind')),.(
    tot_net_generation_mwh = sum(tot_net_generation_mwh),
    tot_fit_net_generation_mwh = sum(tot_fit_net_generation_mwh)), 
    keyby = .(nerc_adj, datetime_utc)
  ]
# Summarizing production by interconnect
ic_gen_dt = 
  elec_gen_dt[!(fuel_category %in% c('solar','wind')),.(
    tot_net_generation_mwh = sum(tot_net_generation_mwh),
    tot_fit_net_generation_mwh = sum(tot_fit_net_generation_mwh)), 
    keyby = .(ic, datetime_utc)
  ]

# First looking at marginal damages 
tot_md_dt = 
  fread(
    here('Data/electricity-generation/output/damage-per-mwh-oge-scc-148.csv')
  ) |>
  melt(
    id.vars = 'plant_id_eia',
    measure.vars = c('md_per_mwh_high','md_per_mwh_low')
  ) |>
  merge(
    plant_info_dt[
      ,.(
      plant_id_eia, 
      nerc_adj,
      stack_height, 
      fuel_category = 
        fcase(
          fuel_category == 'natural_gas', 'Natural Gas',
          fuel_category == 'coal', 'Coal',
          #fuel_category == 'hydro', 'Hydro',
          #fuel_category == 'geothermal', 'Geothermal',
          #fuel_category == 'nuclear', 'Nuclear',
          default = 'Other'
        ) |>
        factor(
          levels = c('Natural Gas','Coal', 'Nuclear','Hydro','Geothermal','Other')
        )
    )],
    by = 'plant_id_eia'
  ) 

damages_hist = 
  ggplot(tot_md_dt) +
  geom_histogram(
    aes(x = value/1000, fill = fuel_category), 
    alpha= 0.5, 
    bins = 100,
    position = 'Identity'
  ) +
  scale_x_log10(labels = label_dollar(accuracy = 0.001)) + #labels = scales::label_dollar()
  theme_minimal() + 
  labs(
    x = 'Marginal Damage ($/KWh)',  
    y = 'Number of Power Plants'
  )+   
  scale_fill_brewer(
    name = 'Fuel Category',
    palette = 'Dark2'
  ) +
  facet_wrap(
    ~ifelse(variable == 'md_per_mwh_high', 'Above median production','Below median production'),
    scales = 'free_y',
    nrow = 2, ncol = 1
  ) + 
  theme(legend.position = 'bottom')
#damages_hist
ggsave(
  plot = damages_hist,
  filename = 'figures/electricity-generation/damages-hist-oge.jpeg',
  width = 6, height = 5,
  bg = 'white'
)  


# Simple predicted vs actual with smoothing -----------------------------------
  ic_gen_dt |> 
  dplyr::group_by(ic)|>
  yardstick::rsq(
    truth = tot_net_generation_mwh, 
    estimate = tot_fit_net_generation_mwh
  )
  # For the interconnect
  model_fit_interconnect_p = 
    ggplot(
      ic_gen_dt, 
      aes(
        x = tot_net_generation_mwh/1e3, 
        y = tot_fit_net_generation_mwh/1e3
      )
    ) + 
    geom_point(alpha = 0.01, color = 'gray30', shape = 19) +
    geom_smooth() +
    geom_abline(intercept = 0, slope = 1, linetype = "dashed") +
    facet_wrap(
      ~ic,
      scales = "free"
    )+ 
    labs(
      x = "Actual Production (GWh)",
      y = "Predicted Production (GWh)"
    ) + 
    theme_minimal() 
  ggsave(
    plot = model_fit_interconnect_p, 
    filename = here('figures/electricity-generation/model-fit-interconnect.jpeg'),
    bg = 'white',
    width = 8, height = 3, units = 'in'
  )
  # ggsave(
  #   plot = model_fit_interconnect_p, 
  #   filename = here('figures/electricity-generation/model-fit-interconnect.pdf'),
  #   device = cairo_pdf,
  #   width = 8, height = 3, units = 'in'
  # )
  model_fit_region_p = 
    ggplot(
      nerc_gen_dt[nerc_adj !='TRE'], 
      aes(
        x = tot_net_generation_mwh/1e3, 
        y = tot_fit_net_generation_mwh/1e3
      )
    ) + 
    geom_point(alpha = 0.01, color = 'gray30', shape = 19) +
    geom_smooth() +
    geom_abline(intercept = 0, slope = 1, linetype = "dashed") +
    facet_wrap(
      ~nerc_adj,
      scales = "free"
    )+ 
    labs(
      x = "Actual Production (GWh)",
      y = "Predicted Production (GWh)"
    ) + 
    theme_minimal() 
  ggsave(
    plot = model_fit_region_p, 
    filename = here('figures/electricity-generation/model-fit-region.jpeg'),
    bg = 'white',
    width = 8, height = 5, units = 'in'
  )
  # ggsave(
  #   plot = model_fit_region_p, 
  #   filename = here('figures/electricity-generation/model-fit-region.pdf'),
  #   device = cairo_pdf,
  #   width = 8, height = 5, units = 'in'
  # )

# Predicted vs actual by hour and season --------------------------------------
  ic_gen_dt[,season:=fcase(
    month(datetime_utc) %in% 3:5, 'Spring',
    month(datetime_utc) %in% 6:8, 'Summer',
    month(datetime_utc) %in% 9:11, 'Fall',
    default = 'Winter'
  )]
  ic_gen_dt[,season:=factor(season, levels = c('Winter','Spring','Summer','Fall'))]
  # Plot of hour of day vs season for interconnection
  model_fit_ic_hour_season_p =
    ic_gen_dt[,.(
      Actual = mean(tot_net_generation_mwh),
      Predicted = mean(tot_fit_net_generation_mwh)), 
      by = .(ic, hour = hour(datetime_utc - hours(5)), season 
    )] |>
    melt(
      #ic_gen_dt,
      id.vars = c('ic','hour','season')
    )|>
    ggplot(aes(x = hour, y = value/1e3, color = variable, linetype = variable)) +
    #geom_smooth(se = FALSE) +
    geom_line(linewidth = 1.1) + 
    facet_grid(
      cols = vars(season), 
      rows = vars(ic), 
      scales = "free"
    ) +
    theme_minimal() + 
    scale_color_brewer(
      palette = "Dark2", 
      labels = c('Actual','Predicted')
    ) + 
    scale_linetype_manual(
      values = c('twodash','solid'), 
      labels = c('Actual','Predicted')
    )+
    labs(
      x = "Hour (EST)", 
      y = "Electricity Production (GWh)",
      linetype = '',color = ''
    ) + 
    theme(legend.position = 'bottom')
  #model_fit_ic_hour_season_p
  ggsave(
    plot = model_fit_ic_hour_season_p,
    filename = here("figures/electricity-generation/model-fit-ic-hour-season.jpeg"),
    bg = 'white',
    width = 8, height = 6
  )
  # ggsave(
  #   plot = model_fit_ic_hour_season_p,
  #   filename = here("figures/electricity-generation/model-fit-ic-hour-season.pdf"),
  #   #device = cairo_pdf,
  #   width = 8, height = 6
  # )
  # Doing the same for regions now 
  nerc_gen_dt[,season:=fcase(
    month(datetime_utc) %in% 3:5, 'Spring',
    month(datetime_utc) %in% 6:8, 'Summer',
    month(datetime_utc) %in% 9:11, 'Fall',
    default = 'Winter'
  )]
  nerc_gen_dt[,season:=factor(season, levels = c('Winter','Spring','Summer','Fall'))]
  model_fit_region_hour_season_p =
    nerc_gen_dt[nerc_adj !='TRE',.(
      Actual = mean(tot_net_generation_mwh),
      Predicted = mean(tot_fit_net_generation_mwh)), 
      by = .(nerc_adj, hour = hour(datetime_utc - hours(5)),season 
      )] |>
    melt(
      id.vars = c('nerc_adj','hour','season')
    )|>
    ggplot(aes(x = hour, y = value/1e3, color = variable, linetype = variable)) +
    geom_line(size = 1.1) + 
    facet_grid(
      cols = vars(season), 
      rows = vars(nerc_adj), 
      scales = "free"
    ) +
    theme_minimal() + 
    scale_color_brewer(
      palette = "Dark2", 
      labels = c('Actual','Predicted')
    ) + 
    scale_linetype_manual(
      values = c('twodash','solid'), 
      labels = c('Actual','Predicted')
    )+
    labs(
      x = "Hour (EST)", 
      y = "Electricity Production (GWh)",
      linetype = '',color = ''
    ) + 
    theme(legend.position = 'bottom')
  #model_fit_region_hour_season_p
  ggsave(
    plot = model_fit_region_hour_season_p,
    filename = here("figures/electricity-generation/model-fit-region-hour-season.jpeg"),
    bg = 'white',
    width = 8, height = 9
  )
#   ggsave(
#     plot = model_fit_region_hour_season_p,
#     filename = here("figures/electricity-generation/model-fit-region-hour-season.pdf"),
# #    device = cairo_pdf,
#     width = 8, height = 9
#   )

# Fuel mix vs load for each interconnection -----------------------------------
  ic_fuel_gen_p = 
    merge(
      elec_gen_dt[,.(
        tot_net_generation_mwh = sum(tot_net_generation_mwh),
        tot_fit_net_generation_mwh = sum(tot_fit_net_generation_mwh)),
        keyby = .(
          ic,
          fuel_category = 
            fcase(
              fuel_category == 'natural_gas', 'Natural Gas',
              fuel_category == 'coal', 'Coal',
              fuel_category == 'hydro', 'Hydro',
              fuel_category == 'nuclear', 'Nuclear',
              #fuel_category == 'geothermal', 'Geothermal',
              default = 'Other'
            ) |>
            factor(levels = c('Natural Gas','Coal', 'Nuclear','Hydro','Other')), 
          datetime_utc
        )
      ],
      nerc_load_dt[,
        .(excess_load = sum(excess_load)), 
        keyby = .(ic, datetime_utc)
      ],
      by = c('ic','datetime_utc')
    )[,':='(
      pct_gen_actual = tot_net_generation_mwh/sum(tot_net_generation_mwh), 
      pct_gen_predicted = tot_fit_net_generation_mwh/sum(tot_fit_net_generation_mwh)), 
      by = .(ic, datetime_utc)
    ] |>
    melt(
      id.vars = c('ic','fuel_category','datetime_utc','excess_load'),
      measure = patterns('pct_gen')
    ) |>
    ggplot(aes(
      x = excess_load/1e3, 
      y = value, 
      linetype = variable,
      color = fuel_category, 
    )) + 
    #geom_point(alpha = 0.2) +
    geom_smooth(se = FALSE) + 
    scale_color_brewer(name = 'Fuel Category',palette = 'Dark2') + 
    scale_y_continuous(labels = scales::label_percent()) +
    labs(
      x = 'Excess Load (GWh)',
      y = 'Percent of Electricity Production'
    ) + 
    scale_linetype_manual(
      name = 'Data',
      labels = c('Actual','Predicted'),
      values = c('twodash','solid')
    ) +
    facet_wrap(~ic, scales = 'free_x')+ 
    guides(
      linetype=guide_legend(keywidth = 4, keyheight = 1, override.aes = list(color = 'black')),
      colour=guide_legend(keywidth = 4, keyheight = 1)
    ) + 
    theme_minimal() + 
    theme(
      legend.position = 'bottom',
      legend.box = 'vertical',
      strip.text.x = element_text(size = 12)
    ) 
  #ic_fuel_gen_p
  ggsave(
    plot = ic_fuel_gen_p,
    filename = here("figures/electricity-generation/model-fit-ic-fuel-mix.jpeg"),
    bg = 'white',
    width = 11, height = 5
  )
  # ggsave(
  #   plot = ic_fuel_gen_p,
  #   filename = here("figures/electricity-generation/model-fit-ic-fuel-mix.pdf"),
  #   #device = cairo_pdf,
  #   width = 11, height = 5
  # )


  # Now fuel mix by region, hopefully explains the variance in 
  region_fuel_gen_p = 
    merge(
      elec_gen_dt[nerc_adj != 'TRE',.(
        tot_net_generation_mwh = sum(tot_net_generation_mwh),
        tot_fit_net_generation_mwh = sum(tot_fit_net_generation_mwh)),
        keyby = .(
          ic,
          nerc_adj,
          fuel_category = 
            fcase(
              fuel_category == 'natural_gas', 'Natural Gas',
              fuel_category == 'coal', 'Coal',
              fuel_category == 'hydro', 'Hydro',
              fuel_category == 'nuclear', 'Nuclear',
              #fuel_category == 'geothermal', 'Geothermal',
              default = 'Other'
            ) |>
            factor(levels = c('Natural Gas','Coal', 'Nuclear','Hydro','Other')), 
          datetime_utc
        )
      ],
      nerc_load_dt[,
        .(excess_load = sum(excess_load)), 
        keyby = .(ic, datetime_utc)
      ],
      by = c('ic','datetime_utc')
    )[,':='(
      pct_gen_actual = tot_net_generation_mwh/sum(tot_net_generation_mwh), 
      pct_gen_predicted = tot_fit_net_generation_mwh/sum(tot_fit_net_generation_mwh)), 
      by = .(ic,nerc_adj, datetime_utc)
    ] |>
    melt(
      id.vars = c('ic','nerc_adj','fuel_category','datetime_utc','excess_load'),
      measure = patterns('pct_gen')
    ) |>
    ggplot(aes(
      x = excess_load/1e3, 
      y = value, 
      linetype = variable,
      color = fuel_category, 
    )) + 
    #geom_point(alpha = 0.2) +
    geom_smooth(se = FALSE) + 
    scale_color_brewer(name = 'Fuel Category',palette = 'Dark2') + 
    scale_y_continuous(labels = scales::label_percent()) +
    labs(
      x = 'Excess Load (GWh)',
      y = 'Percent of Electricity Production'
    ) + 
    scale_linetype_manual(
      name = 'Data',
      labels = c('Actual','Predicted'),
      values = c('twodash','solid')
    ) +
    facet_wrap(~nerc_adj, scales = 'free_x')+ 
    guides(
      linetype=guide_legend(keywidth = 4, keyheight = 1, override.aes = list(color = 'black')),
      colour=guide_legend(keywidth = 4, keyheight = 1)
    ) + 
    theme_minimal() + 
    theme(
      legend.position = 'bottom',
      legend.box = 'vertical'
    ) 
  #region_fuel_gen_p
  ggsave(
    plot = region_fuel_gen_p,
    filename = here("figures/electricity-generation/model-fit-region-fuel-mix.jpeg"),
    bg = 'white',
    width = 11, height = 8
  )
  # ggsave(
  #   plot = region_fuel_gen_p,
  #   filename = here("figures/electricity-generation/model-fit-region-fuel-mix.pdf"),
  #   #device = cairo_pdf,
  #   width = 11, height = 8
  # )


# Plant level fit plots -------------------------------------------------------

summarize_plant_gen_dt_plant = function(nerc_adj_in){
  # Getting file paths for all plants in region
  nerc_plant_gen_fp = 
    plant_gen_fp[
      str_extract(plant_gen_fp, '(?<=-)\\d*(?=\\.fst)') |> as.integer()
      %in% plant_info_dt[nerc_adj == nerc_adj_in]$plant_id_eia
    ]
  # Summarizing generation by time, ic, nerc, fuel cat
  elec_gen_dt = 
    map_dfr(
      nerc_plant_gen_fp, 
      \(fp_in){
        read.fst(
          path = fp_in,
          as.data.table = TRUE
        ) |>
        merge(
          plant_info_dt, 
          by = 'plant_id_eia'
        )
      }
    )[,.(
      tot_net_generation_mwh = mean(net_generation_mwh),
      tot_fit_net_generation_mwh = mean(fit_net_generation_mwh),
      tot_prod_hours = mean(net_generation_mwh > 0),
      tot_fit_prod_hours = mean(net_generation_mwh > 0)
      ), 
      keyby = .(
        plant_id_eia,
        ic = factor(fcase(
            nerc_adj %in% c('CAL','WECC'), 'West',
            nerc_adj == 'TRE', 'Texas',
            default = 'East'
          ), levels = c('West','Texas','East')),
        nerc_adj, 
        fuel_category
      )
    ]
      
  return(elec_gen_dt)
}

plant_gen_dt = 
  map_dfr(
    c('CAL','MRO','NPCC','RFC','SERC','TRE','WECC'),
    summarize_plant_gen_dt_plant
  )

# Total production hours predicted vs actual 
plant_hours_p =
    plant_gen_dt[,.(
      Actual = tot_prod_hours,
      Predicted = tot_fit_prod_hours,
      tot_mwh = tot_net_generation_mwh), 
      by = .(plant_id_eia, ic)
    ] |>
    ggplot(aes(x = Actual, y = Predicted, size = tot_mwh, weight = tot_mwh)) + 
    geom_point(alpha = 0.3, color = 'gray30', shape = 19) +
    geom_smooth() +
    geom_abline(intercept = 0, slope = 1, linetype = "dashed") +
    # facet_wrap(
    #   ~ic,
    #   scales = "free"
    # )+ 
    scale_x_continuous(
      name = 'Actual Production Hours', 
      labels = scales::label_percent()
    )+ 
    scale_y_continuous(
      name = 'Predicted Production Hours', 
      labels = scales::label_percent()
    ) + 
    theme_minimal() +
    guides(size = 'none')
 # plant_hours_p
  ggsave(
    plot = plant_hours_p,
    filename = here("figures/electricity-generation/model-fit-plant-hour.jpeg"),
    bg = 'white',
    width = 8, height = 6
  )
  # ggsave(
  #   plot = plant_hours_p,
  #   filename = here("figures/electricity-generation/model-fit-plant-hour.pdf"),
  #   device = cairo_pdf,
  #   width = 8, height = 6
  # )
  
  # Total production by plant ---------
  plant_prod_p =
    plant_gen_dt[,.(
      Actual = tot_net_generation_mwh,
      Predicted = tot_fit_net_generation_mwh), 
      by = .(plant_id_eia, ic)
    ] |>
    ggplot(aes(x = Actual, y = Predicted)) + 
    geom_point(alpha = 0.3) + #alpha = 0.3, color = 'gray30', shape = 19
    #geom_smooth() +
    geom_abline(intercept = 0, slope = 1, linetype = "dashed") +
    # facet_wrap(
    #   ~ic,
    #   scales = "free"
    # )+ 
    scale_x_continuous(
      name = 'Actual Production (MWh)' ,
      labels = scales::label_number(scale = 1, big.mark = ',')
    )+ 
    scale_y_continuous(
      name = 'Predicted Production (MWh)',
      labels = scales::label_number(scale = 1, big.mark = ',')
    ) + 
    theme_minimal()
  #plant_prod_p
  ggsave(
    plot = plant_prod_p,
    filename = here("figures/electricity-generation/model-fit-plant-prod.jpeg"),
    bg = 'white',
    width = 8, height = 6
  )
  # ggsave(
  #   plot = plant_prod_p,
  #   filename = here("figures/electricity-generation/model-fit-plant-prod.pdf"),
  #   device = cairo_pdf,
  #   width = 8, height = 6
  # )
    
  # Total production by plant-month ---------
  summarize_plant_gen_dt_plant_month = function(nerc_adj_in){
  # Getting file paths for all plants in region
  nerc_plant_gen_fp = 
    plant_gen_fp[
      str_extract(plant_gen_fp, '(?<=-)\\d*(?=\\.fst)') |> as.integer()
      %in% plant_info_dt[nerc_adj == nerc_adj_in]$plant_id_eia
    ]
  # Summarizing generation by time, ic, nerc, fuel cat
  elec_gen_dt = 
    map_dfr(
      nerc_plant_gen_fp, 
      \(fp_in){
        read.fst(
          path = fp_in,
          as.data.table = TRUE
        ) |>
        merge(
          plant_info_dt, 
          by = 'plant_id_eia'
        )
      }
    )[,.(
      tot_net_generation_mwh = mean(net_generation_mwh),
      tot_fit_net_generation_mwh = mean(fit_net_generation_mwh),
      tot_prod_hours = mean(net_generation_mwh > 0),
      tot_fit_prod_hours = mean(net_generation_mwh > 0)
      ), 
      keyby = .(
        plant_id_eia,
        ic = factor(fcase(
            nerc_adj %in% c('CAL','WECC'), 'West',
            nerc_adj == 'TRE', 'Texas',
            default = 'East'
          ), levels = c('West','Texas','East')),
        nerc_adj, 
        fuel_category, 
        month = month(datetime_utc)
      )
    ]
      
  return(elec_gen_dt)
}

plant_month_gen_dt = 
  map_dfr(
    c('CAL','MRO','NPCC','RFC','SERC','TRE','WECC'),
    summarize_plant_gen_dt_plant_month
  )


  plant_month_prod_p =
    plant_month_gen_dt[,.(
      Actual = tot_net_generation_mwh,
      Predicted = tot_fit_net_generation_mwh), 
      by = .(plant_id_eia, month, ic)
    ] |>
    ggplot(aes(x = Actual, y = Predicted)) + 
    geom_point(alpha = 0.3) + #alpha = 0.3, color = 'gray30', shape = 19
    #geom_smooth() +
    geom_abline(intercept = 0, slope = 1, linetype = "dashed") +
    #facet_wrap(
    #  ~ic,
    #  scales = "free"
    #)+ 
    scale_x_continuous(
      name = 'Actual Production (MWh)' ,
      labels = scales::label_number(scale = 1, big.mark = ',')
    )+ 
    scale_y_continuous(
      name = 'Predicted Production (MWh)',
      labels = scales::label_number(scale = 1, big.mark = ',')
    ) + 
    theme_minimal()
  #plant_month_prod_p
  ggsave(
    plot = plant_month_prod_p,
    filename = here("figures/electricity-generation/model-fit-plant-prod-month.jpeg"),
    bg = 'white',
    width = 8, height = 6
  )
  # ggsave(
  #   plot = plant_month_prod_p,
  #   filename = here("figures/electricity-generation/model-fit-plant-prod-month.pdf"),
  # #  device = cairo_pdf,
  #   width = 8, height = 6
  # )

